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Abstract 



We study the long-time behavior of stochastic models with an ab- 
sorbing state, conditioned on survival. For a large class of processes, 
in which saturation prevents unlimited growth, statistical properties 
of the surviving sample attain time-independent limiting values. We 
may then define a quasi- stationary probability distribution as one in 
which the ratios Pn{t)/Pmit) (for any pair of nonabsorbing states n and 
m), are time- independent. This is not a true stationary distribution, 
since the overall normalization decays as probability flows irreversibly 
to the absorbing state. We construct quasi-stationary solutions for the 
contact process on a complete graph, the Malthus-Verhulst process, 
Schlogl's second model, and the voter model on a complete graph. We 
also construct the master equation and quasi-stationary state in a two- 
site approximation for the contact process, and for a pair of competing 
Malthus-Verhulst processes. 
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I. INTRODUCTION 



Stochastic processes with an absorbing state arise frequently in statistical physics 
epidemiology 0, and related fields. On one hand we have autocatalytic pro- 
cesses (in lasers, chemical reactions, or surface catalysis, for example), and popula- 
tion models (in genetic or epidemiological modeling), in which a population of n > 
individuals goes permanently extinct when the absorbing state, n = 0, is reached. 
On the other hand, phase transitions to an absorbing state in spatially extended 
systems , exemplified by the contact process 0], are currently of great interest in 
connection with self-organized criticality the transition to turbulence 0, and, 
more generally, issues of universality in nonequilibrium critical phenomena 
Such models are fequently studied using deterministic mean-field equations, Monte 
Carlo simulation, and renormalization group analyses. 

While fiuctuations play an essential role in the vicinity of the absorbing state, 
exact solution of the master equation is not, in general, feasible. A useful tool in this 
situation is van Kampen's fi-expansion in inverse powers of the system size, about 
the deterministic, macroscopic solution. It is clearly desirable to develop additional 
methods for analyzing processes with absorbing states. In this work we study models 
whose macroscopic equation admits a nonabsorbing (or active) stationary state, but 
which, for a finite system size, must always end up in the absorbing state. We show 
that the quasi- stationary distribution for such a system is readily constructed, and 
that it provides a wealth of information about its behavior. In fact, simulations of 
"stationary" properties of lattice models with an absorbing state actually study the 
quasi- stationary regime, since the only true stationary state (for a finite system) is 
the absorbing one. 

For models without spatial structure, such as uniformly distributed populations, 
well stirred chemical reactors, or networks in which each node communicates equally 
with all others, our results provide a complete description of long-time properties, 
conditioned on survival. In models with spatial structure, typified by nearest- 
neighbor interactions on a lattice, as in the contact process, the description in terms 
of one or a few variables corresponds to a mean-field theory that cannot adequately 
capture critical fiuctuations. Our approach does, nonetheless, allow one to put some 
of the fiuctations "back into" mean-field theory, and so we are able to study finite 
size effects and moment ratios that are beyond the grasp of simpler approaches. For 
one-step models involving a single variable, the computational demands of our ap- 
proach are trivial; more complicated processes can be analyzed using straightforward 
numerical procedures. 

The idea of the quasi-stationary state or distribution is quite simple. Consider 
a Markov chain > in continuous time, with n = absorbing, and such that the 
macroscopic or rate-equation description (i.e., neglecting fluctuations) admits an 
active stationary state [ngt > 0). Then in many cases one expects the probability 
distribution Pn{t) to bifurcate, after some initial transient, into two components. 



2 



one with all weight on the absorbing state, the other, g„(t), concentrated near the 
macroscopic value rist, such that p„ f» for n ~ 0. The survival probability (for the 
process not to have fallen into the absorbing state), is P(t) = J2n>i Qnit). We study 
the case in which g„ has attained a time-independent form, in which its only time 
dependence is through the overall factor Ps{t). Since the probability distribution is 
supposed to bifurcate, we may write 

Pnit) = Q{t)5nfl + qn{t) (1) 

where Q{t) = \ — P(t), and by definition, go(^) = 0. The quasi- stationary distribu- 
tion (QS) is then defined via the condition: 

p„(t) = P(t)p„, {n > 1), (2) 

where the p„ are time independent. (Again, pQ = 0.) Since P(t) is the survival 
probability we have the normalization: 

EPn = l- (3) 

n>l 

Clearly, not every process with an absorbing state admits a QS distribution. The 
exponential decay process (with transition rates Wn,m = ^Sn,m-i) and the birth- 
and-death process {Wn,m = m6n,m-i + ^^Sn,m+i) , for example, do not. While we 
will not try to determine rigorously the conditions under which a QS distribution 
exists, it seems reasonable to suppose that the macroscopic equation for the process 
in question admits an active stationary solution. All of the examples considered 
satisfy this condition. The voter model represents a special case, in which the 
macroscopic equation provides no information on the fate of the process, and the 
QS distribution bears no resemblance to the stationary (absorbing) distribution. 

The balance of this paper is organized as follows. In Sec. 2 we illustrate the eval- 
uation of the QS distribution, and the analysis of associated statistical properties, 
for the simple case of the contact process on a complete graph. Sec. 3 is concerned 
with the closely related Malthus-Verhulst process (MVP). Both the contract process 
and the MVP exhibit a continuous phase transition in the infinite-size limit. In Sees. 
4 and 5 we study a models with a discontinuous phase transition: the second Schlogl 
process, and the voter model. In Sec. 6 we return to the contact process, this time 
in a two-site approximation, and in Sec. 7 study a pair of competing MVPs. Sec. 8 
contains a summary and discussion of our results. 



II. CONTACT PROCESS ON A COMPLETE GRAPH 

In the contact process (CP) @J^, each site of a lattice is either occupied {cTiit) = 
1), or vacant (cr,;(t) = 0). Transitions from cTj = 1 to cTj = occur at a rate of 
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unity, independent of the neighboring sites. The reverse can only occur if at least 
one neighbor is occupied: the transition from cTj = to cTj = 1 occurs at a rate 
of Am, where m is the fraction of nearest neighbors of site i that are occupied; 
thus the state cxj = for all i is absorbing. The stationary order parameter p (the 
fraction of occupied sites), is zero for A < Ac. It is easy to show that in mean- field 
approximation (i.e., treating occupancy of sites as statistically independent events, 
and supposing that the density p{t) is spatially uniform), the density follows the 
rate equation: 

|^(A-l).-A/. (4) 

This equation predicts a continuous phase transition (from p = Otop = l — A^^ in 
the stationary state) at Ac = 1. While this is qualitatively correct, Ac is in general 
larger than 1; Ac = 3.29785(2) in one dimension, for example. 

A mean-field description such as Eq. (1) does not, of course, yield correct critical 
exponents for dimension d < = 4, but that is not of concern here. Instead, we 
would like to study the simplest stochastic process for which Eq. (||) represents the 
macroscopic limit, (i.e., the limit of an infinite system, in which p{t) is a deterministic 
variable). By studying such a process we can recover the fluctuations ignored in the 
macroscopic equation. Since the process has an absorbing state, its limiting (t — oo) 
distribution is trivial {pn(t) — Sn,o), but we can study the long-time properties, 
conditioned on survival, using the QS distribution. 

The stochastic birth-and-death process whose macroscopic limit is given by Eq. 
(1) is the contact process on a complete graph, in which the rate for a vacant site to 
turn occupied is A times the fraction of all sites that are occupied, rather than the 
fraction of nearest neighbors. Since each site interacts equally with all others, all 
pairs of sites are neighbors, defining a complete graph. Given its connections with 
epidemic modeling and percolation, this model has also been called "percolitis" 0. 

The state of the process is specified by a single variable n: the number of occupied 
or infected sites. This is a one-step process with nonzero transition rates 

Wn-l,n = n (5) 
71 

W^n+i,n = A-(fi-n) (6) 

on a graph of Q sites. These expressions yield the master equation 
dpn 



dt 



{n + l)pn+i + {n-l)X{n-n + l)pn-i- 1 + X{^1 - n) np„ , (7) 



where A = X/Q. Since P{t) = J2n>iPnit), we have dP/dt = — Pi(t) in the contact 
process, and under the QS hypothesis, 



4 



If we insert Eq. (Q) in the master equation, then, dividing through by P{t) we 
obtain 



X{n-l)[n-n+ + + - [X{n-n) + l]np^+p,p^ = 0, (9) 

for n > 1. For a given n, Eq. (^) furnishes a relation for p^+i in terms oi p^ and 
Letting Un = n[X{Q — n) + 1], we have 

P2 = -Pi) (10) 



and 



Pn = - 

n 



r 



{Un-1 - Pi)Pn-l + (^-2)(l^-2)Ap„_2 



'11^ 



for n = 2, fi. Thus the QS distribution is completely determined once Pi is known; 
the latter is determined via the normalization condition, Eq. (j^). 

In practice the following iterative scheme converges very quickly. Starting with 
a guess for p^ (unity, for instance), one uses the recurrence relations to find the 
corresponding and then evaluates S = J2n=iPn- The procedure is then repeated 
using p[ = Pi/S; after a modest number of steps, S converges to unity, at which 
point the p„ represent the QS distribution. 

We have confirmed that the master equation for the contact process on a complete 
graph does in fact attain the QS distribution at long times. We integrate the master 



equation using a fourth-order Runge-Kutta scheme 0], and stop the integration 
when the mean population conditioned on survival, {n)s = J2nWn(t)/P(t), changes 
by less than 10~^° per time step. The resulting distribution, = pn(t)/P(t) is 
compared with the QS distribution furnished via the recurrence relations in figure 
1; they are identical. 

Having constructed the QS distribution and verified that it indeed represents the 
long-time distribution (conditioned on survival), as given by the master equation, we 
now examine some properties of the QS state. Figure 2 shows the quasi-stationary 
population density p = n/Q versus A for various system sizes. This plot looks 
strikingly similar to finite-size plots of the order parameter at a continuous phase 
transition. In particular, it is evident that the A-dependence is smooth for any finite 
system, but that on increasing Q, the function p(A, Q) approaches a singular limit. 

On finite-dimensional lattices the value of the order parameter at the critical 
point is expected to scale as p(Ac, f^) ~ where P and z/_|_ are the critical 



exponents governing the order parameter and the divergence of the correlation length 
0. Here we find p{Xc,^) ~ 0.685(^-^/2 ^^^^^ results for fl = 500 - 10^). We note 
that the mean- field exponent values (3=1 and = 1/2 for the contact process 
would lead one to expect p{\c, ^) ~ The reason for the difference would seem 
to be that on a complete graph (where each site is, in effect, a unit distance from all 
others), the correlation length is undefined, and so finite-size scaling ideas do not 
apply. Instead, we can understand the exponent -1/2 by noting that the microscopic 
equation (^, gives p = at the critical point, so that n is purely a fiuctuation. The 
central limit theorem then requires n ~ fi^/^. Away from the critical point, A = 1, 
the finite-size correction to the mean-field solution, p = 1 — A^^, is 0{Q~^). 

Another property of interest at a continuous phase transition is the moment ratio 
m = (?7,2)/(?T,)2. This quantity is analogous to Binder's reduced fourth cumulant 
||1 1|| , at an equilibrium critical point: the curves m{X,Q) for various Q cross near 
Ac (the crossings approach Ac), so that m assumes a universal value rric at the 
critical point. For the basic contact process (and other models in the universality 
class of directed percolation), rric = 1.174 and 1.326 in one and two dimensions, 
respectively |T^. In figure 3 we plot m for the quasi-stationary contact process on 
a complete graph. Qualitatively, the behavior is similar to that observed on finite- 
dimensional lattices, with m(A) becoming ever steeper with increasing system size; 
\dm/d\\x^ ~ f2^/2 for large ^l. The crossings approach Ac, as expected. We find that 
m(l, n) ~ 1.660 + C((]-i/2) ig^^gg ^ 

As noted above, the survival probability decays as dP{t)/dt = —pi{t). We may 
therefore interpret as the quasi-stationary decay rate. The QS lifetime, r = l/pi, 
grows as exp [const. Q{X — Ac)] for A > Ac, as can be seen from figure 4, where we 
plot In r versus A for several system sizes. At the critical point, r ~ fi^'^. 

It may appear surprising that we are able to discuss quasi-stationary properties 
for A < Ac, where an active stationary state does not exist, even in the thermo- 
dynamic limit. The existence of nontrivial QS properties for A < Ac is in fact a 
finite-size effect, analogous to a nonzero magnetization in the Ising model above 
the critical temperature, on a finite lattice. For A > Ac the QS state converges, as 
Q oo, to the true stationary state, while for A < Ac it converges to the absorbing 
state. But since the properties of any finite system are nonsingular, we must expect 
a smooth decay of, for example, the density, as A — >• 0. 



A. Pseudo-stationary distribution 

The contact process on a finite graph (be it a regular lattice, a complete graph, 
or small-world network) has only the absorbing state as its true stationary state. 
It is nonetheless instructive to examine the result of forcing a time-independent 
solution on the master equation. This is easily done by introducing the generating 
function F{z,t) = J2nPn(t)z"'. The master equation, Eq.(|^), is equivalent to the 
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partial differential equation, 



dF 



- BF ^ d ( dF\ 



(12) 



Setting dF/dt = and letting G{z) = dF/dz, a first integration yields 

G{z) = Cz^'-^e'/^' , (13) 

where C is a constant. Integrating once more gives 

oo Q-n 

Fps{z) = CJ:y^ r- (14) 

While the corresponding distribution, 

Pps,m = TTT Tj (15) 

(ii — m)!m 

yields dpm/dt = when inserted in the master equation, it is not a valid solution 
since it cannot be normalized, and does not satisfy the correct boundary condition 
at m = 0. For this reason we call Pps,m. a pseudo-stationary distribution. Similar 
pseudo-stationary distributions have been derived by Munoz from the Fokker-Planck 
equation for continuous versions of the contact process and of a model with multi- 



plicative noise [T^ 



If we treat Pps,m as a proper probabihty distribution, restricting it to m = 1, f2, 
and normalizing it on this set, then the resulting distribution turns out to be very 
close to the quasi-stationary one when is negligible in the vicinity of m = 0. 
Figure 5 (for Q = 100) shows that the two distributions are essentially identical 
for A = 2 (as they are for larger A values), and for A = 1.5, except for small n. 
For A = 1.2, however, the distributions are radically different; in particular, we see 
that the pseudo-stationary distribution does not respect the absorbing boundary at 
n = 0. 



In summary, the pseudo-stationary distribution provides a simple and useful 
approximation to the true QS distribution for A (and Q) sufficiently large that 

~ for m ~ 1. (In this case the survival probability decays extremely slowly, 
rendering the hypothesis of a stationary distribution - while strictly incorrect - at 
least reasonable.) As we approach the critical value, however, the region near n = 
assumes dominance, and Pps bears little relation to the quasi-stationary distribution. 
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III. MALTHUS-VERHULST PROCESS 



The Malthus-Verhulst process (MVP) is a birth-and-death type process, n{t), 
very similar to the contact process on a complete graph, but with the difference 
that n may assume any nonnegative integer value. The nonzero transition rates for 
this one-step process are: 

Wn-i,n = fJ'n + ^n{n-l) , (16) 

Wn+i,n = Xn . (17) 

In the first expression, Q represents the system size, as discussed by van Kampen [|l|. 
In the MVP there is no lattice structure, and no fixed limit to the population size, 
as there is in the contact process; rather, growth is limited by the term oc z/, which 
represents competition between individuals for access to finite resources. (Thus 
u represents an intrinsic competition parameter, and (z//f2)n(n — 1) is the rate of 
interactions between pairs of individuals in a system of size Q.) The saturation 
effect that permits a nontrivial quasi- stationary state is imposed in the death term, 
whereas it appears in the birth term, in the contact process. 



Letting n = Qp, the macroscopic equation for the MVP is 

{X-p)p-up' (18) 



dt 

which has the stationary solution p = (A — /i)/z/. The master equation is 
dpn 



dt 



[p + un]{n+l)pn+i + {n-l)Xpn-i - [X + p + z/(n-l)]np„ , (19) 



with = I'/Q. Inserting the quasistationary form of the probability distribution, 
one readily obtains the recurrence relation 

- _ (gn-l - /iPl)Pn^l - - 2)p„_2 , . 

n[p + z/(r2 — 1)J 

where 

qn = [\ + p + V{n-l)]n (21) 

The QS distribution can be found using the same iterative procedure as for the 
contact process. (We cut off the distribution at a certain M such that p<10^^° for 
n> M.) 
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The quasi-stationary properties of the MVP are quite similar to those of the 
contact process on a complete graph. Figure 6 shows the dependence of p on A for 
system sizes Q = 100, 10^ and 10^. (We fix /i = z/ = 1 in all numerical studies of 
the MVP.) The density grows linearly with A for A > Ac = /i, unlike in the contact 
process, where the macroscopic saturation term is also proportional to the birth 
rate, A. As in the contact process on a complete graph, the lifetime r ~ fi^/^ at Ac, 
and the QS density decays ~ fi"-*^/^. 

The moment ratios m (figure 6, inset) show the same qualitative trend as in 
the CP. The asymptotic value of m at the critical point, moreover, is the same as 
in the CP: we find m ~ 1.660 — 0.814fi~^/^ for large il. This identity of m values 
suggests that the limiting QS probability distributions at Ac are the same for the two 
processes. This is verified numerically; for Q = 10^, for example, the distributions 
differ by less than about 2 x 10^^. 

In fact, the critical QS probability distribution for the MVP (and by extension 
for the contact process) enjoys the scaling property 

p„ = -l^f{n/Vn) . (22) 



The scaling function / may be found using a method that parallels van Kampen's 
analysis of the master equation |jl|. First, we treat n in Eq. (|20|) as a continuous 
variable, expanding Pn±i to second order. Setting A = /x at the critical point, we 
obtain 

- [(2/i + V)n + Vv?] — — + [2yU + z7n(l + n)] — + 2Vnp = -fiPiP . (23) 
/ dn^ an 



Next we make the change of variables y = n/n^^'^, and f{y) = Vl^^'^p^. This yields 
an expansion in inverse powers of fi^/^; at lowest order (f2~^) we find 

^,y^ + [2/. + uy']^ + 2uyf = -/x/(0)/(y) . (24) 



The initial value /(O) must be chosen such that fdy = 1. 

To integrate Eq. (p^) numerically we let g{y) = df/dy, and note that g{0) = 
— [/(0)]^/2, and dg/dy\y=o = — (2z///i)/(0). Setting p = = 1, numerical integration 
yields the scaling function shown in figure 7, which agrees well with the distributions 
for the MVP and the contact process on a complete graph, in the limit of large system 
size. [Note that the change of variables used to derive Eq. (^) yields the very same 
equation when applied to the recursion relation for the contact process, Eq. (|^), at 
its critical point, A = 1.] The asymptotic scaling function f{y) yields the moment 
ratio m = 1.660063, again in agreement with results from the recurrence relations. 
(Note that m = /(O)/ (y), as may be verified directly from Eq. (0).) 
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Despite some superficial differences, the MVP and the contact process on a 
complete graph show the same type of phase transition (continuous, with mean- 
field-like critical exponents), in the limit Q — * oo, and the same scaling properties 
at their respective critical points. 

As in the case of the contact process, one may define a pseudo-stationary distri- 
bution by demanding that the r.h.s. of the master equation be zero. The resulting 
expression, 

Pps,n = C— I ^ ) -TT (25) 

is very similar to that of the contact process on a complete graph, and again repre- 
sents a good approximation to the quasi-stationary distribution if ~ for n ~ 1. 



IV. SCHLOGL'S SECOND MODEL 



The processes studied in the preceding sections exhibit a continuous phase tran- 
sition; next we consider a process with a discontinuous transition. Schlogl's second 
model [|14| describes the set of autocatalytic reactions A — 0, 2 A SA, and 

—i> 2/1 in a well stirred system. Since the growth process is quadratic in the den- 
sity (rather than linear, as in the contact and MV processes), a low-density active 
state is unstable, and we expect a discontinuous transition to a state with a nonzero 
density as the growth rate is increased. The nonvanishing transition rates are: 



V 



n— l.n 



-l)(n-2) 



(26) 



+l,n 



A 



-n(n — \\ 



(27) 



These yield the macroscopic equation: 

dp 



dt 



HP + \p — up 



(2^ 



which admits an active stationary solution 

1 
2 



Ps 



A + JA2 - Apu 



(29) 



for A > 2^//Iz7 (The stationary density jumps from zero to y/pi' at the transition. 
It is known that in a spatially extended system the transition is actually continu- 
ous for d < A |jl6|,|15[, but here we treat the well stirred system, which exhibits a 
discontinuous phase transition.) 
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Since this is a one-step process, it is a simple matter to derive the recurrence 
relations for the quasi-stationary probabilities: 



n[fi + v{n-l){n-2)] ^ ' 

where 

qn=[fi + J{n-l)+V{n-l){n-2)]n, (31) 

with X = X/Q and V = v/VL^. These may be solved using the same method as for 
the contact process or the MVP. 

In the following numerical example, we fix /i = z/ = 1, and study the neighbor- 
hood of the transition at A = 2. Figure 8 shows the quasi-stationary density near the 
transition; with increasing VL the density approaches the discontinuous stationary 
solution, Eq. p9|). The density goes from 0{1/VL) to 0(1), in a transition region of 
width ~ The moment ratio m, at the transition, appears to approach unity 

(from above) slowly, as ^ oo; data for < 5 x lO'' yield m ~ 1 + 2.64fi°-^^^. 
Similarly, the density at the transition approaches its limiting value from below: 
p ~ 1 — 1.5f2°'^^^. The lifetime at the transition appears to grows as r ~ 
Essentially the same power laws (but with different prefactors) are found at the 
transition for fi = 2 and u = 1/2. 

The probability distribution is bimodal in the neighborhood of the transition, a 
hallmark of a discontinuous transition. This is evident in figure 9, where the relative 
amplitudes of the two peaks, one at n = 1, the other near n ^ 1000, shift rapidly 
near the transition. The bimodal distribution of course presages a hysteresis loop 
when Q oo. 



V. VOTER MODEL ON A COMPLETE GRAPH 



A somewhat simpler model exhibiting a discontinuous phase transition is the 
voter model [|17|,|18|. It differs from the systems considered previously in that it 
possesses two absorbing states. In the voter model each site on a lattice is in one 
of two states, A or B; each A — B nearest-neighbor pair has a unit rate to change 
to the state AA (with probability uj) or to BB (with probability 1 — uo). Let n{t) 
denote the number of sites in state A; both n = and n = Vt are absorbing, on 
a lattice of VL sites. Starting from a random initial condition with equal densities 
of A and B, the stationary probability distribution switches discontinuously from 
5nfl to 5nn as u is increased through 1/2. Of particular interest is the coarsening 
dynamics for uj = 1/2 [Q. In this case {n) is constant, and the probabilities of the 
two absorbing states are determined by n{{]). 
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In this brief section we study the model on a complete graph, for the case uj = 
1/2. The master equation is 



This yields the macroscopic equation dp/dt = 0, reflecting the dominance of fluctu- 
ations in this process. From the terms for n = and n = Q we find that in the QS 
state, the survival probability obeys 

^^ = -^(l-^-')(Pi+Pn), (33) 



allowing us to write the following relation for the QS distribution: 

(34) 



The solution is the uniform distribution, p^ = — 1) for n = 1, ...,Q — 1 (here, 
Po = Pn = 0). 

In the case of the voter model, the QS distribution looks very different from the 
stationary distribution, even in the infinite-fi limit. This is because the process has 
no active stationary state: on a complete graph it must always fluctuate into one of 
the absorbing states. If we exclude all trials that hit one of the boundaries, we are 
in effect looking at a random walker confined to an interval, whose distribution is 
asymptotically uniform. 



VI. TWO-SITE APPROXIMATION FOR THE CONTACT PROCESS 



In this section we return to the contact process, constructing a stochastic descrip- 
tion in terms of nearest-neighbor pairs of sites, on a ring of Q sites. The idea of the 
approximation is similar to that of two-site or pair mean-field dynamic mean-field 
theory [0, except that in the present case we treat the number of occupied sites, 
n, and the number of doubly-occupied nearest-neighbor pairs, z, as stochastic vari- 
ables. When necessary, we invoke the usual factorization of three-site probabilities 
in terms of those for two sites, in order to construct the transition rates. 

We begin by establishing the range of allowed values for z. Using '1' and '0' to 
represent occupied and vacant sites, respectively, z is the number of (11) nearest- 
neighbor (NN) pairs. Let w represent the number of (10) NN pairs. (For obvious 
reasons, the number of (01) pairs is again w.) Note that w is not an independent 
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variable, since each 1 is followed wither by a or another 1, yielding n = z + w. 
Similarly, the number of (00) pairs v is given hj v = Q — 2n + z. (Here we used 
n + 2w + z = Q.) Evidently v is restricted to the set {0, Q} while w must be in 
{0, Q/2} (We assume fl even.) These considerations imply certain limits for z on 
a ring of Q sites, as listed in table 1. 



Table 1. Allowed values for z in the CP on a ring. 



n 


z 


0, 1 





2,...,n/2 


0,...,n-l 




2n — Q, ...,n — 1 


Q 


Q 



Next we construct the transition rates Wn',z';n,z- From a given state {n, z) there 
are at most five possible transitions, viz., to the states {n', z') = {n—1, z), {n—1, z—1), 
{n—l,z — 2), {n + 1, z + 1), and {n + 1, z + 2). The first three represent annihilation 
of a particle with, resp. zero, one or two occupied nearest neighbors, while the last 
two represent creation of a particle at a site with one, or two occupied neighbors, 
resp. The rate of the transition (n, z) {n—1, z) is the number of isolated particles, 
that is, triplets of the form (010). This number is not determined completely by n 
and z] we estimate it as the number of (01) pairs times the conditional probability 
of a vacant site given an occupied neighbor: w x [w/n). This estimate, however, 
is obviously wrong when z = n — 1 and n > 2, in which case there are no (010) 
triplets. Thus we have 



n—l,z:n,z 



0, 



z = n — 1, n > 1 
otherwise 



(35) 



By similar arguments, we find: 



Wn-l,z-l;n,z 

Wn-l,z-2;n,z 

,z+l;n,z 
,z+2;n,z 



[n — z 



0, 



z = n — 1, n > 1 
I — , otherwise 

Z = l 



n — 2, z = n — 1, n = 3, — 1 
— , otherwise 



A, 



z 



n — 1 



X{n — z) ^^'^^ , otherwise 
0, 



z = n — 1 
otherwise 



(36) 

(37) 

(38) 
(39) 



In the macroscopic limit, the master equation defined by these rates yields a pair 
of coupled equations for the particle density p and the pair density ( = z/Q 0. The 
critical value Ac = 2 in this approximation. 
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Since this is not a one-step process, we cannot obtain the QS distribution via 
recurrence relations, as was done in the previous examples. The most obvious ap- 
proach - direct integration of the master equation - works reasonably well for modest 
system sizes {Q < 200 or so), but becomes very slow for large Q. This is due in part 
to the large number of equations (~ and to the fact that the transition rates 
are 0{fl), forcing one to use a time step ~ fi^^, to avoid numerical instability. A 
much more efficient iterative scheme that converges directly to the QS distribution 



is described in Ref. |]T9[. The basic idea of this approach is that in the QS regime, 

Pn = , (40) 

Wn - ro 



for any nonabsorbing state n, where r„ = J2n' Wn^n'Pn'{t) is the flux of probability 
into state n {tq is the flux into the absorbing state), and Wn = J2n'Wn',n is the 
total transition rate exiting n. The iterative scheme starts with an arbitrary 
distribution p„ normalized on n > 1 (or, in the multivariate case, on the set of 
nonabsorbing states). Next one evaluates 

T 

Pn = aPn + {l-a) — , (41) 



where a is a parameter between zero and one. The new distribution p'^ is then 
normalized, and inserted in the r.h.s. of Eq. (^Tj), and the process repeated until it 
converges. For large Q the time to reach the QS distribution is three or more orders 
of magnitude smaller than that required for numerical integration. (In the present 
case convergence is enhanced if one uses a small value for a, e.g., 0.01. A further 
economy is realized by noting that, near Ac, Pnz is extremely small for n ~ 
allowing one to truncate the distribution. For Q = 1280 and A = 2, for example, 
~ 10-1^ for n = 350.) 

We studied the QS properties of the contact process on a ring of Q = 10, 20, 
40,..., 2560 sites, using the pair approximation. Figure 10 shows that the dependence 
of the particle density and the moment ratio on A and on system size is qualitatively 
similar to that found in the complete graph case. The scaling laws for the density 
(oc ^7^^/^) and the lifetime (oc fi^/^) at the critical point are the same as for the CP 
on a complete graph. The critical values of the moment ratio m at Ac, extrapolated 
to infinite system size, yield rric = 1.653, nearly the same as for the CP on a complete 
graph. The critical probability distribution p:^ again shows a scaling collapse (figure 
11), but the scaling function differs somewhat from that for the complete graph. 
Figure 12 shows the joint QS distribution ^ at the critical point, for Q = 40. For 
each n there is a most probable number, z*{n), of pairs; z* grows linearly with n. 
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VII. COMPETING MALTHUS-VERHULST PROCESSES 



As a final example we consider a population consisting of two types, A and 
B, that evolve according to the MVP analyzed in Sec. 3. (The two subpopulations 
may, for example, possess different alleles of a certain gene, in a species with asexual 
reproduction.) A and B interact via two mechanisms: (1) each individual competes 
with all others (regardless of type) for access to resources; (2) on reproducing, A 
mutates to B with probability q, and vice- versa. Thus we have the following set of 
transition rates: 

Wn^-l,nB;nA,nB = I^AUa + ^^^(^^+^^-1) , (42) 



WnA,nB-l;nA,nB = t^BTlB + ^^£(^^ + ^5-1) , (43) 



WnA+l,nB;nA,nB = ^APTlA + As^ris , (44) 
WnA,nB+l;nA,nB = >^BpnB + A^^n^ , (45) 



where p = 1 — q. 

Evidently the state UA — nB — O is absorbing. For mutation probability q — 0, 
riA = is also absorbing (similarly ub — 0), and in the QS state only one type 
is present, i.e., PnAUB = if both ua and ub are nonzero. This prompts us to 
investigate the effect of a nonzero mutation rate: Will the predominance of one 
species persist when mutation is possible? To address this issue we study the order 
parameter 

A ^ (46) 

{ua + Ub) 

in the QS state. Figure 13 shows how A varies with g, for various system sizes. The 
competing MVPs are both at the critical point, Xi — jii. (For simphcity we take all 
rates, Aj, /Xj, and z/, equal to unity). The order parameter decreases with increasing 
mutation probability, as expected, and with increasing system size. From analysis 
of A at various g-values, as a function of f), it appears that lima^oo A(g, VL] Ac) = 
for any nonzero q. In other words, an arbitrarily small mutation rate restores the 
symmetry of the process. We observe qualitatively similar behavior in A(g, Q) above 
the critical point (A > 1). At the critical point, the decay of A(g, Q) with Vl is 
slower than exponential, but faster than a power law; it can be fit approximately by 
a stretched exponential ~ exp[— c where the factor c depends on q (see figure 
13 inset). 
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VIII. DISCUSSION 



We study quasi- stationary probability distributions for Markov processes exhib- 
ing a phase transition between an active and an absorbing state. The QS distri- 
bution, which represents the long-time limit of the probability distribution, condi- 
tioned on survival, yields a complete statistical description of the process in this 
limit. While the QS distribution can always, in principle, be found by integrating 
the master equation, we present much more efficient methods for its generation: 
recurrence relations for one-step processes, and an iterative scheme for multi-step 
or multivariate processes. 

Some insight into the nature of the quasi-stationary distribution is afforded by 
the considering the equation (in the notation of Sec. VI), for n > 1, 

^ = -WnQn + rn + roQn , (47) 



where Vn = J2n' Wn,n'(ln' IS the flux of probability into state n > 0. (As usual we 
take state as absorbing.) Without the final term on the right-hand side, this is 
simply the master equation. Including this term, however, the equation has as its 
stationary solution the QS distribution, In other words, the QS distribution 
corresponds to the stationary state of a "process" in which all fiux to the absorbing 
state is immediately reinserted into the non-absorbing subspace. The portion alloted 
to state n is equal to its probability, g„. 

Thus the QS distribution does not result from imposing a simple boundary con- 
dition (such as periodic, or refiection at the origin), on the original master equation. 
Since Eq. (0) is nonlinear, it is not a master equation. It may happen, nevertheless, 
that the QS distribution for a certain process (possessing an absorbing state), is also 
the stationary distribution for some other process, as in fact was found for the voter 
model. 

In this work we study QS distributions for processes with an absorbing state, 
but without spatial structure. We find that the Malthus-Verhulst process and the 
contact process on a complete graph have the same scaling properties near their 
respective critical points. A study of Schlogl's second model illustruates how a dis- 
continuous phase transition emerges in the infinite-size limit, associated with a rapid 
change in the QS distribution, which is bimodal in the vicinity of the transition. The 
voter model on a complete graph serves as a negative example: it has no active sta- 
tionary state, and the QS distribution, which is uniform on the set of allowed density 
values, is very different from the true long-time distribution, even in the infinite-size 
limit. Our two final examples are bivariate processes. The pair approximation to 
the contact process shows scaling properties that are very similar to those of the 
complete-graph version. A study of a pair of competing Malthus-Verhulst processes 
reveals that mutation causes a radical change in the QS distribution. 
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Several directions for future work on quasi-stationary distributions can be men- 
tioned. First, it should be possible to derive continuum equations, analogous to the 
Fokker-Planck equation, to describe the QS distribution in the limit of large system 
size. Second, analysis of QS properties for a series of system sizes promises to be a 
useful method for studying lattice models possessing an absorbing state. Through 
finite-size scaling analyses and extrapolation procedures, it should be possible to 
extract useful estimates of critical properties for such models. Finally, numerical 
study of more elaborate models used in epidemiology and population dynamics may 
yield a better understanding of such processes. 
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FIGURE CAPTIONS 



Figure 1. Probability distributions for the contact process on a complete graph, 
L — 100 sites. The curves (connecting values at integer points) represent the long- 
time limit of the master equation, normahzed to the survival probability P{t). Points 
represent the quasi-stationary distribution determined by the recurrence relations. 
Left curve: A = 1.2; right: A = 2. 

Figure 2. Quasi-stationary population density versus A for the contact process on 
a complete graph. System sizes 100, 200, 500, 1000, and 2000, top to bottom on the 
left-hand side. The inset is a detail of the critical region, including the curve for the 
infinite-size limit, p — 1 — A~^. 

Figure 3. Quasi-stationary moment ratio m versus A for the contact process on a 
complete graph. System sizes 20, 50, 100, 200, 500, and 1000, in order of increasing 
steepness. 

Figure 4. Quasi-stationary lifetime r versus A for the contact process on a complete 
graph. System sizes 100, 500, 1000, and 5000, upper to lower on left-hand side. 

Figure 5. Quasi-stationary distribution (solid line) and pseudo-stationary distri- 
bution (points) for the contact process on a complete graph, L = 100 sites for A = 2 
(upper), 1.5 (middle), and 1.2 (lower). The inset in the middle panel is a detail of 
the region near n — 0. 

Figure 6. Quasi-stationary density versus A for the Malthus-Verhulst process. 
System sizes 100, 10^, and 10^, upper to lower on the left-hand side. The inset 
shows moment ratio m versus A for the same system sizes. 

Figure 7. Scaling plot of the quasi-stationary distribution, / = f^^^^Pn versus 
y — n/fl^^"^, at the critical point: +: MVP, Vl = 10^; o: contact process on complete 
graph, fl = 10^; solid curve: asymptotic scaling function f{y) obtained via numerical 
integration. 

Figure 8. Quasi-stationary density versus A for the second Schlogl process with 
fjL — u — l. System sizes 100, 500, 1000, and 5000, (in order of increasing steepness). 
The dashed curve is the mean-field solution, which is discontinuous at A = 2. 

Figure 9. Quasi-stationary distribution for the second Schlogl process with p — 
v — 1, system size Q, — 1000. A = 1.98, 2, and 2.02 (upper to lower on left-hand 
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side) . 

Figure 10. Quasi-stationary density versus A for the contact process in the pair 
approximation. System sizes 20, 40, 80, 160, 320, and 640. The inset shows moment 
ratio m versus A for the same system sizes. 

Figure 11. ScaUng plot of the quasi- stationary distribution: P* = f^^^^p^ versus 
n* = at the critical point of the contact process. Open circles: pair approx- 

imation, Q = 640; sohd hue: pair approximation, Q = 2560; filled circles: contact 
process on complete graph, Q = 1000. 

Figure 12. Quasi-stationary probabihty distribution for the critical contact process 
in the pair approximation, system size D, — AO. 

Figure 13. Quasi-stationary order parameter A versus mutation rate q for a pair of 
competing MVPs at the critical point. System sizes = 20, 50, 100, 200, 500, 1000 
(upper to lower). The inset shows In A versus 0°'^ for (upper to lower) q — 0.02, 
0.05, and 0.1. 
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